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CONVERSION  FACTORS,  U.S.  CUSTOMARY  TO  METRIC  (SI)  UNITS  OF  MEASUREMENT 

U.S.  customary  units  of  measureme.t  used  in  this  report  can  be  converted  to 
metric  (SI)  units  as  follows: 


Multiply 

by 

To  obtain 

inches 

25.4 

millimeters 

2.54 

centimeters 

square  inches 

6.452 

square  centimeters 

cubic  inches 

16.39 

cubic  centimeters 

feet 

30.48 

centimeters 

0.3048 

meters 

square  feet 

0.0929 

square  meters 

cubic  feet 

0.0283 

cubic  meters 

yards 

0.9144 

meters 

square  yards 

0.836 

square  meters 

cubic  yards 

0.7646 

cubic  meters 

miles 

1.6093 

kilometers 

square  miles 

259.0 

hectares 

knots 

1.852 

kilometers  per  hour 

acres 

0.4047 

hectares 

foot-pounds 

1.3558 

newton  meters 

millibars 

1.0197  x  10-3 

kilograms  per  square  centimeter 

ounces 

28.35 

grams 

pounds 

453.6 

grams 

0,4536 

kilograms 

ton,  long 

1.0160 

metric  tons 

ton,  short 

0.9072 

metric  tons 

degrees  (angle) 

0.01745 

radians 

Fahrenheit  degrees 

5/9 

Celsius  degrees  or  Kelvins1 

lTo  obtain  Celsius 

(C) 

temperature  readings 

from  Fahrenheit  (F)  readings. 

use  formula:  C  *  (5/9) 

(F  -32). 

To  obtain  Kelvin 

(K) 

readings,  use  formula: 

K  -  (5/9)  (F  -32)  +  273.15. 
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SYMBOLS  AND  DEFINITIONS 


ai  ,bl 


12 

d 

db 

E 

F 

fn 

GB,GBP 

g 

Hb 

id 

k 

L 

i 
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N 

n 

Pts 

P1 

P 

Q12 

R 

S12 


XHF 

v 

z 

6 

Y 

0 

Ad 

Af 

At 

(tf 


Fourier  series  coefficients 

distance  from  bottom  to  pressure  sensors 

wave  celerity 

cospectrum  value 

total  water  depth 

breaking  wave  depth 

wave  energy  density 

complex  Fourier  coefficient 

discrete  frequency  value 

ratio  of  rms  breaking  wave  height  to  breaking  wave  depth 
acceleration  of  gravity 
breaking  wave  height 
counting  indexes 

dynamic  pressure  response  factor 
wave  number 
wavelength 
sensor  spacing 

Index  to  account  for  gage  number 
total  number  of  discrete  data  points 

frequency  nimiber ,  argument  of  Fourier  series  coefficients 

longshore  energy  flux  at  the  surf  line 

pressure  time-series  values 

dynamic  pressure 

quad-spectrum  value 

ratio  of  unwlndowmd  energy  density  to  windowed  energy  density 

complex  cross-spectnm  value 

length  of  time  series  record 

high  frequency  cutoff  period 

weighting  coefricient 

water  surface  elevation 

gage  orientation  angle 

specific  weight  of  seawater 

wave  direction 

average  mean  depth  of  water  overlaying  pressure  sensors 
frequency  step 
time  step 

angular  wave  frequency 


COMPUTER  ALGORITHM  TO  CALCULATE  LONGSHORE  ENERGY  FLUX  AND 
WAVE  DIRECTION  FROM  A  TWO  PRESSURE  SENSOR  ARRAY 


by 

Todd  L.  Walton,  Jr.  and  Robert  G.  Dean 

I.  INTRODUCTION 

The  documented  (FORTRAN  IV  programing  language)  computer  program  discussed 
in  this  report  was  originally  written  as  part  of  the  Coastal  Engineering 
Research  Center’s  (CERC)  Longshore  Sand  Transport  Research  Program  and  was  used 
in  analysis  of  wave  data  collected  at  Channel  Islands  Harbor  in  conjunction 
with  a  study  of  sand  transport  at  Channel  Islands  Harbor  as  discussed  in  Bruno, 
et  al.  (1981). 

The  program  performs  the  basic  analysis  of  two  wave  gage  pressure  records 
necessary  to  compute  wave  direction  and  wave  energy  at  a  given  frequency  and 
computes  the  longshore  energy  flux  used  in  sand  transport  for  the  entire  energy 
spectrum  of  the  wave  record.  This  program  uses  linear  wave  theory  for  the  wave 
transformation  process  and  includes  the  assumption  of  straight  and  parallel 
bottom  contours  necessary  for  application  of  Snell's  law  of  refraction. 

Necessary  steps  in  the  analysis  of  the  wave  data  are  presented  in  Sections 
II  and  III  of  this  report.  Subroutines  are  discussed  and  sample  outputs  for 
some  wave  records  from  the  Channel  Islands  wave  gage  pressure  sensor  pair  are 
given. 

The  program  presently  accepts  data  in  the  standard  CERC  magnetic-tape  for¬ 
mat  where  record  lengths  consist  of  4,100  values.  The  ffrst  four  values  are 
the  gage  number  and  the  date-time  group,  and  the  remaining  4,096  values  are 
the  pressures  recorded  in  thousandths  of  a  foot  (head)  of  water  at  0. 25-second 
intervals.  Should  other  input  data  be  available,  the  program  could  easily 
be  modified  to  accept  the  data  by  simple  changes  in  the  main  program  and  in 
subroutines  BUF  and  SWITCH. 

Sample  outputs  have  been  presented  for  real  wave  data;  some  wave  direc¬ 
tional  information  cannot  be  obtained  for  all  frequencies  because  the  spectral 
information  at  some  frequencies  is  ill-conditioned.  Hie  percent  of  energy  for 
which  this  problem  occurs  is  a  small  part  of  the  energy  (usually  <3  percent)  of 
the  entire  spectrum  and  is  insignificant  in  energy-flux  computations.  Reasons 
for  this  feature  are  discussed  later. 

II.  METHODOLOGY 

Calculating  the  longshore  energy  flux  at  breaking  required  the  following 
steps : 

(1)  Calculation  of  the  frequency -by-frequency  wave  direction  and 
energy  at  the  location  of  the  wave  gages; 

(2)  determination  of  the  breaking  wave  depth; 

(3)  transformation  of  the  wave  spectrum  to  the  "breaker"  line, 
including  shoaling  and  refraction  effects;  and 

(4)  computation  of  “W  the  longshore  energy  flux  at  the 
surfline. 

Each  of  the  steps  is  described  below. 


1 .  Calculation  of  Wave  Direction  and  Energy  Spectrum  at  Wave  Gages, 

As  noted  previously,  each  of  the  input  time-series  pressure  records  con¬ 
sists  of  4,096  data  points  with  a  time  increment  of  0.25  second.  To  reduce 
computational  costs,  modified  time  series  are  formed  for  analysis  by  averaging 
four  adjacent  data  points.  These  new  time  series  contain  1,024  data  points 
spaced  at  1.0-second  intervals.  This  increases  the  aliasing  period  from  0.5 
to  2.0  seconds;  however,  this  is  justified  as  the  pressure  response  factor 
for  a  water  depth  of  6  meters  and  a  wave  period  of  2  seconds  is  approximately 
0.005. 

The  time  series  are  analyzed  using  a  standard  fast  Fourier  transform  (FFT) 
program  to  determine  the  coefficients.  For  example,  for  pressure  time  series 
from  gage  1 


Pl(j>  -  I  l®i(n)  -  ib!(n)]  exp  U) 

n*  0  V  / 

in  which  i  -  /-l  and  N  is  the  total  number  of  data  points,  T/At  *  1,024, 
where  T  is  the  time  series  record  length  of  1,024  seconds,  At  the  time 
increment  of  1  second  between  samples,  and  j  a  discrete  time  tj  where  tj  * 
discrete  time  value  »  jAt.  The  FFT  coefficients  are  defined  in  terms  of  the 
pressure  time  series  as 


ai (n)  -  ib^n) 


N-l 

I  l 

N  4  n 

j-0 


(2) 


where  the  argument  "n"  of  the  Fourier  coefficients  a(n)  and  b(n)  speci¬ 
fies  the  quantity  to  be  a  discrete  function  of  wave  frequency,  fR,  where 
fn,  a  discrete  frequency  value,  is  nAf  (where  Af  ■  1/T)  and  the  a^O)  term 
represents  the  mean  value  of  the  time-series  pressure  record  for  wave  gage  1. 
Similar  relationships  exist  for  wave  gage  2.  In  calculating  the  FFT  coeffi¬ 
cients,  there  are  several  options  that  may  be  employed  in  an  attempt  to  reduce 
spectral  leakage  which  arises  due  to  representing  an  aperiodic  time  series  by 
a  periodic  series.  A  large  number  of  possible  data  windows  (weighting  func¬ 
tions  for  data)  have  been  developed  to  reduce  the  adverse  effects  of  spectral 
leakage  (Harris,  1974).  These  can  be  expressed  in  the  form  of  a  weighting 
function  w( j ) ,  such  that  the  modified  time  series  p’(j)  is  of  the  form 

p’(j)  -  w(j)  p(j) 

in  which  p(j)  is  the  digitized  measured  pressure  value  at  time  tj  «  jAt,  and 
w(j)  a  weighting  function.  A  characteristic  of  these  weighting  functions  is 
that  they  are  equal  to  unity  at  the  midpoint  of  the  time  series  and  decrease 
to  a  lesser  value  near  the  two  ends.  In  the  present  program,  a  " cosine  bell** 
weighting  function  is  used;  however,  through  comparisons  of  P »8  with  and 
without  this  function,  it  was  established  that  the  effect  of  the  weighting 
function  was  minimal  (<5  percent).  The  cosine  bell  weighting  function  is 
expressed  by 


w(j) 


1.0  -  cos 


(3) 
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It  is  clear  that  the  application  of  a  weighting  function  will  reduce  the  total 
energy  in  the  record.  This  effect  is  partly  compensated  for  by  the  following 
equation: 


P”(j) 


■J<E- 

V  <p,2> 


p'(j) 


(4) 


thereby  ensuring  the  same  total  energy  in  the  altered  and  original  time 
series,  where  <p2>  is  the  mean  square  value  of  the  original  time  series 
and  <p,2>  the  mean  square  value  of  the  weighted  time  series.  It  is  the 
altered  time  series  p"(j)  that  is  subjected  to  FFT  analysis.  The  primes 
will  be  dropped  hereafter  for  convenience.  The  average  mean  depth  of  water 
overlying  the  pressure  sensors,  Ad,  is  obtained  by  averaging  the  m  time 
series  to  obtain  am(0).  For  two  separate  time  series  records,  m  *  1,  2  (wave 
gages  1  and  2), 


Ad  =  0.5  [a  L(0)  +  a2( 0)  ]  (5) 

The  total  water  depth,  d,  is  the  sum  of  Ad  and  the  distance,  B,  o^ 
the  pressure  sensors  above  the  bottom  (in  later  examples  B  =  0.76  meter). 

Each  FFT  pressure  coefficient  is  transformed  to  a  water  surface  displace¬ 
ment  coefficient  by  the  following  linear  wave  theory  relationship  discussed  in 
the  Shore  Protection  Manual  (SPM)  (see  Ch.  2,  U.S.  Army,  Corps  of  Engineers, 
Coastal  Engineering  Research  Center,  1977): 

water  surface  coefficients  dynamic  pressure  coefficients 

1 

[afn),  b  (n) ]  -  -  [a ra(n),  bm(n)l  (6) 

in  which  the  subscripts  n  and  p  denote  water  surface  and  dynamic  pressure 
coefficients,  respectively.  The  factor 


Kz(n) 


cosh  k(n)  B 
cosh  k(n)  d 


(7) 


where  y  is  the  specific  weight  of  fluid  (seawater)  and  is  included  when 
pressure  coefficients  are  in  normal  units  of  pressure  (i.e.,  N/M2  or  equiva¬ 
lent).  In  equation  (7),  B  represents  the  distance  of  the  pressure  sensors 
above  the  bottom  and  k(n)  is  the  wave  number  associated  with  the  angular 
frequency,  u>(n)  -  (2wnAf),  as  obtained  from  the  linear  wave  theory  dispersion 
relationship 


u>(n)2  «  gk(n)  tanh  k(n)  d  (8) 

One  of  the  disadvantages  of  measuring  waves  with  near-bottom  pressure  sen¬ 
sors  is  evident  by  examining  equations  (6)  and  (7).  For  the  higher  frequen¬ 
cies  (shorter  wave  periods)  K^n)  is  very  small  which  means  that  the  higher 
frequency  waves  result  in  very  small  pressure  fluctuations  near  the  sea  floor. 
Thus,  to  avoid  contaminating  the  calculated  water  surface  displacements,  it  is 


usually  necessary  to  apply  a  high  frequency  cutoff,  above  which  the  pressure 
contributions  are  discarded.  The  proper  selection  of  this  high  frequency 
cutoff  depends  on  the  signal  to  noise  characteristics  of  the  pressure  sensor 
and  the  signal  conditioning  system.  In  the  present  program,  the  high  fre¬ 
quency  cutoff  was  established  at  a  wave  period  of  3.0  seconds.  Wave  gage 
analyses  by  Thompson  (1980)  have  shown  that  a  3.0-second  high  frequency 
spectral  cutoff  value  provides  reasonable  estimates  of  total  wave  energy  at 
west  coast  (U.S.)  locations. 

Denoting  hereafter  the  FFT  coefficients  for  the  water  surface  as  a(n) 
and  b(n) ,  it  is  noted  that  the  coefficients  have  the  following  properties: 

N-x 


<n2> 

-  1  [a2(n)  + 

b2(n) J 

(9) 

and 

n*i 

•(* 

+  n) 

(10) 

‘(f 

+  n) 

i- -*(*-■) 

(11) 

and  thus 

N/  2 

<n2> 

=2  l  (a2(n) 

+  b2  (n)] 

(12) 

n-l 


Thus,  the  total  (kinetic  and  potential)  energy  E(n)  associated  with  a  par¬ 
ticular  wave  frequency  component,  n,  is 

E(n)  =  2y[a2(n)  +  b2(n) ]  (13) 

Now  consider  two  wave  or  pressure  sensors  located  at  (x^ ,  y  )  and  (x2  ,  y^) 
(see  Fig.  1).  The  results  will  be  developed  considering  discrete  frequencies. 


Figure  1.  Definition  sketch  for  two  sensor  array. 


10 


The  water  surface  displacement  consistent  with  the  assumption  of  one 
direction  per  frequency  is 


N-l 

n(x,  y,  j)  -  l  F(n)  exp  Ulnu^t  -  k^n)  x  -  ky(n)  y)} 

n-o 


N-l 

*  I 


n-o 


[a(n)  -  ib(n)]  exp 


(14) 


where  is  the  primary  analysis  frequency  (*  2ir/record  length  ■  2w/T  - 

2irAf),  and  0(n)  the  direction  of  wave  propagation  at  frequency  w(n)  -  no^  • 
The  wave  number  components,  k^n)  and  ky(n),  are  expressed  in  terms  of  the 
wave  number,  k(n),  and  wave  direction,  0(n),  as 


kx(n)  *  k(n)  cos  0(n) 

(15) 

ky(n)  -  k(n)  sin  0(n) 

(16) 

The  cross  spectrum,  S12(n),  of  the  two  measured  water  surface  displacements 
(or  dynamic  pressures)  is  given  by 


S12(n)  «  | F(n) | 2  {exp  -  i  [k(n)  cos  0(n)(x2  -  xx) 
+  k(n)  sin  0(n)(y2  -  y^]} 


(17) 


Denoting  the  separation  distance  and  angle  as  i  and  0,  respectively,  the 
cross  spectrum  can  be  expressed  as  (see  Fig.  1) 

Sj2(n)  *  |F(n)|2  {cos  [k(n)  icos  (0(n)  -  0) ]  -  i  sin  [k(n)  icos  (0(n)  -  0)]} 

*  cospectrun  (n)  -  i  quad-spectrum  (n)  (18) 

*  Ci2(n)  -  iQi2(n) 


Thus,  from  equation  (18),  the  wave  direction  0(n)  associated  with  each  wave 
frequency  can  be  expressed  as 


(19) 


The  above  relationship  has  two  roots,  one  of  which  must  be  selected  based 
on  physical  considerations  of  the  most  likely  direction  of  wave  propagation. 
In  the  present  case,  assuming  no  wave  reflection  from  the  beach,  the  ambiguity 
in  wave  direction  is  ruled  out;  for  wave  sensors  nearly  parallel  to  the  beach, 
the  minus  sign  in  equation  (19)  is  appropriate. 
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There  are  two  conditions  for  which  it  was  not  possible  to  calculate  the 
wave  directions  0(n) .  These  include  poorly  conditioned  wave  data,  presumably 
due  to  spectral  leakage,  and  spatial  aliasing  due  to  large  separation  distance 
between  the  two  gages.  If  the  data  are  poorly  conditioned  for  determining 
wave  direction,  the  absolute  value  of  the  quantity  within  the  brackets  {-}  in 
equation  ( 19)  may  exceed  unity,  a  physically  impossible  condition  since  the 
extreme  values  of  the  cosine  function  are  ±1*  This  tends  to  occur  for  the 
extremely  long  waves  for  which  the  energy  is  small  and  the  value  of  k(n)  is 
also  small,  the  latter  tending  to  result  in  large  values  of  the  bracketed 
quantity.  The  percentage  of  energy  for  which  this  .cndition  occurred  in  the 
analysis  of  one  year’s  wave  data  collected  at  Channel  Islands  Harbor  was 
relatively  small,  averaging  2  to  3  percent  with  a  maximum  of  approximately 
10  percent.  The  second  condition  is  related  to  spatial  aliasing  and  requires 
that  one-half  the  wavelength  be  equal  to  or  greater  than  the  projection  of  the 
wave  gage  separation  distance  in  the  direction  of  wave  propagation.  Referring 
to  Figure  1 , 


L  >  21  {cos ( 0(n)  -  ej Jmax  (20) 

which  indicates  that  for  the  least  adverse  effects  of  spatial  aliasing,  the 
gages  should  be  on  an  alinement  parallel  to  the  dominant  orientation  of  the 
wave  crests.  As  will  be  discussed  later,  in  calculating  P^g  an  attempt  was 
made  to  account  for  this  effect  of  aliasing  by  augmenting  the  calculated 
values,  illustrated  as  follows  by 


ET0T 

(p£s^cm  ”  (P£s^c 


(21) 


in  which  the  subscripts  c  and  cm  indicate  calculated  and  calculated  modi¬ 
fied,  respectively.  et0T  and  E  represent  the  total  wave  energy  values 
and  the  wave  energy  not  affected  by  spatial  aliasing  or  poorly  conditioned 
wave  data,  respectively.  The  total  wave  energy  is  that  energy  in  the  wave 
spectrum  below  the  high  frequency  spectral  cutoff  value. 

2.  Transformation  of  Wave  Spectrum  to  Breaker  Line. 

At  this  stage,  the  wave  energy  and  wave  direction  in  the  vicinity  of  the 
gages  are  determined.  These  values  are  then  transformed  to  the  breaker  line 
accounting  for  wave  refraction  and  shoaling. 

To  determine  the  wave  breaking  depth,  the  onshore-directed  energy  flux  is 
calculated  in  accordance  with  the  expression  (based  on  Snell’s  law  of  refrac¬ 
tion)  and  equated  to  an  equivalent  expressed  in  terms  of  wave  characteristics 
at  breaking . 


N/  2 

Onshore  energy  flux  *  £  y2  [a(n)2  +  b(n)2]  Cg(n)  cos  0(n) 

n-i 


cgb  cos 


(22) 
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*•  ■ 


PIIWJ.IU| 


Assuming  that  the  breaking  wave  angle,  Ob,  is  small,  that  the  waves  will 
break  under  shallow-water  conditions,  and  that  the  ratio  of  breaking  wave 
height  to  depth  is  a  constant,  the  breaking  wave  height,  H^,  is  then  given 
by 


(N/2  „  )  /< 

2B  \°  * 2 

|  l  16  [a(n)2  +  b(n)2]  Cg(n)  cos  0(n)j0,14^- 

i) 

(23) 

where  GB  is  the  ratio  of  root-mean-square  (rms)  breaking  wave  height  to 
breaking  depth,  GB  =  H^/d^  (here  assumed  GB  -  0.78)*  With  the  breaking  depth 
known,  each  wave  component  is  transformed  to  shore  accounting  for  both  wave 
refraction  and  shoaling  based  on  linear  wave  theory. 

Wave  refraction  is  in  accordance  with  Snell's  law  and  the  assumption 
that  straight  and  parallel  contours  existed  between  the  gage  and  breaking 
locations 


0b(n) 


*  sin”1 


Cu 

— 2 —  sin  0  (n) 
Cr(n)J 


(24) 


where  C  is  linear  wave  celerity  (see  the  SPM,  Ch.  2)  in  which  the  r  sub¬ 
scripts  denote  the  "reference  (gage)'*  location. 

With  the  wave  energy  and  direction  now  known  at  the  breaker  line,  the 
value  of  the  longshore  energy  flux,  (Pj£S)cm>  is  readily  determined 


<Pis>cm  “ 


is'c 


(  N/2 

R  (2y  l  [a2(n)  +  b2(n)]b  [Cg(n)]b[cos  0(n)  sin  0(n)J 


(25) 


n-l 


in  which  the  factor  R  is  given  by  the  ratio 


ET0T 

R  =  - 

E 


as  defined  in  and  discussed  in  relation  to  equation  (21). 


III.  MAIN  PROGRAM  DOCUMENTATION 

The  detailed  programing  steps  in  analysis  for  the  longshore  energy  flux, 
(P^g)cmi  (which  in  this  program  is  calculated  in  terms  of  rms  wave  height)  are 
presented  in  this  section.  Program  steps  are  numbered  to  correspond  to  areas 
in  the  program  listing  where  computations  are  carried  out.  A  program  listing 
with  corresponding  numbered  steps  follows  the  program  documentation.  Note 
that  preceding  text  has  used  the  indexes  j  and  n  for  time  and  frequency, 
respectively,  while  the  program  which  follows  uses  the  index  I  for  both  time 
and  frequency.  A  listing  of  the  main  program  is  presented  in  Figure  2. 
Program  steps  are  as  follows  and  refer  to  numbered  parts  of  main  program 
listing : 
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10 


15 


FHUCKAH  *P|CT UNPUT . OUTPUT* TAPt5»lNPUT,TARtSPOUTPUT# MPEG) 
t  COMPUTOR  ALGORITHM  10  CALCULATE  LONORMORt  ENERGY  FLU*  FACTOR  AND  WAVE 

C  DIRECTION  FOR  TWO  P*t*»y«*  SENSOR  ARRAY 

C  HAIM  PROGRAM 

C  PROGRAM  II  PRESENTLY  lit  UP  TO  TAKE  A  U*t  SERIES  OF  111*  FGlMf  IN  **1* 
OlHfMliON  C ( S 1 2* 

DIMENSION  FIR(  |0|i)  *F  tl(  IPS*}  «F|R(  1 0l«)*F21Moa«> 

U I Hfc M* JON  IICMA(llt) •FH00lQ($l*)tTHETA(5ia) 

DIMENSION  Ct2(51i)«0ll(f l*) 

OIHCNljON  M(t02«) 

DIMENSION  CC(5l2)  •IINTHSCMII 
HlAL  H(AM| t Ht AN2 
LOGICAL  1*0 
DATA  END/. FALSE./ 

101  f  UMMATU0C2**FS.*)) 


ao 


a5 


10 


15 


AO 


AS 


10 


55 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

L 

c 

t 

c 

c 

c 

c 

t 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


OtUNmON$*Fl*€D  VARIABLES 

A»t*FONENTlAL  POWER  DEFINING  NUMBER  OF  TIME  SEMES  PO INTM  <  ***K  ) 

SMPACING  BETWEEN  WAVE  GAGES  (FEET) 

OtlTTwTlHC  STf  P  BETWEEN  POINTS  In  AVERAGED  TIME  SCRieS  (IECONOS) 

BETA  ■ANGLE  DIFFERENCE  BETWEEN  WAVE  CAGE  ALIGNMENT  ANp  SHORE LINE (RADI ANS) 
»LOPE«SLOPE  OF  BEACH  AT  POINT  OF  WAVE  BREAKINq 
GAMMAbsPECIFIC  WEIGHT  OF  FLUID  (LBS/FTMl) 

BbUISTaNCE  OF  PRESSURE  SENSORS  ABOVE  BOTTOM  (FEU) 

(^ACCELERATION  OF  GRAVITY  (FEET/Sec***) 

C0»RATIO  BREAKING  WAVE  HUGHT/UEPTM  FOR  LIN( AR  THEORY  COMPUTATION 

GBPbNaTIO  BREAKING  WAVE  HEIGHT/OEPTH  FOR  LINEAR  THEORY  COMPUTATION  OF 
•utFR  DEPTH  CIYEN  BREAKING  WAVE  HEIGHT 


DEFINITION® •FLOAT  I  NO  VARIABLES 
AVGl>AvERAG e  OF  time  SERIESI 
AVG2* AVERAGE  OF  TIME  SERIES  2 
C(1)WWAVI  celerity 
ctauucosPccTRA  or  &ERiisi*a 
CbaBREAKlNC  wave  CELERITY 
CG(I)«GROUP  wave  CfiURlTV 

CNTLU)««096  POINT  TIME  SERIES  BEFORE  AVERAGING 

0|PTH«dIPTH  OF  waT|R  AT  GAGE  SUL  FROm  AVERAGES  OF  GAGES  )  AND  2 

F\1(I)»VJN0EF  XNED/CQWPLEX  IMAGINARY  PORTION  OF  TRANSFORM 

FtR(l)ST|*E  SERIES  DATA  GAGE I/COHPLE*  REAL  PORTION  Of  TRANSFORM 

F2i(I)*UNDEFtNeQ/CQMPLE*  IMAGINARY  PO«TlON  OF  TRANSFORM 

F2N ( I  )«T JMf  SCRIES  DATA  GAGE2/C0MPLE*  RE*l  PORTION  OF  TRANSFORM 

FMU0SU( UPTIME  SERIES  AMPLITUDE  moouLOS  SQUARED 

HBBBKEAKING  wave  height 

1A(I)B5000  POINT  DATA  CROUP  ANU  TIME  SERIES  RECORD 

plneg«negative  ccnthibutiun  to  longshore  energy  flu*  factor 
PLNETaNET  LONGSHORE  tNERGT  FLU*  FACTOR 

PLP0S*P0SITI YE  CUN T H 1  HU T ION  TO  LONGSHORE  ENERGY  FLU*  FACTOR 
U1*(I)«<JUaD9PECTWA  OF  SERIES  1*2 

r«scaling  factor  for  scaling  up  energy  of  nonusable 

GohTionB  OF  DIRECTIONAL  IPtCIRA 
MAIIUURATIO  OF  EN|HUY/WINUQWEU  ENERGY  fur  gage  I 
rat io2«r a t xu  of  cnergy/winoo^eo  energy  fur  gage  2 

RE AL( I )■ 1020  POINT  TIME  SERIES  AFTER  AVERAGING 


*0 


♦5 


TO 


C  RNbRaTIO  of  group  wave  cileritt  to  WAVE  CELERITY 

C  RShFRUbPERCENT  OF  ENERGY  BEYUNU  SPACIAL  ALIASING  FREQUENCY 

C  HSLFRUb PERCENT  OF  ENERGY  BELOW  LOW  FREguENCv  CUTOFF 

C  HSUQUaPfMCCNT  OF  INCOHERENT  ENERGY 

C  SHFREUPSUH  UF  ENERGY  WITH  FREQUENCIES  ABOVE  SPACIAL 
C  ALIAS INO  FREQUENCY  CUTOFF 

C  SHG2PSUMMATIUN  OF  ONSHORE  ENERGY  FLUX 

C  SIGMA(2)«RA01AL  FREQUENCY 

C  SLFHEUbSUM  OF  ENIRQY  wjTh  FREQUENCIES  BELOW  LOW  F REGuENCY  CUTOFF 
C  SUOORSUM  OF  ENERGY  with  FREQUENCIES  HAVING  INCOHERENT  RAVE  OIRECTION 
C  SUH|N«SUH  UF  ENERGY 

C  SUmimum  OF  SQUARES  UF  TIME  SERIES  t  WITHOUT  AVERAGE 

C  SUHt'SUM  OF  SQUARES  OF  TIME  SERIES  2  WITHOUT  AVERAGE 

C  SUWF1MUW  OF  SQUARES  OF  TIRE  SERIES  1  w|TM  AVEHAOI 

C  SUMF2PSUH  OF  SQUARES  UF  TIME  SERIES  2  WITH  AVERAGE 

C  TPWAVE  PERIOD 

t  TMETA(I)RWAV|  direction  in  RaOIANS 


Figure  2 •  Listing  of  main  program. 
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Ti 


§0 


a* 


to 


C  T Kfc T ABaBHE A* INQ  *AvE  ANQLl 

C  n8uhi«ium  or  squares  of  om  nindon  modified  time  hmii  i 

C  OUMfaiUM  OF  SQUARES  OF  DM*  *1NQ0*  MUDIFUO  TI«I  IfDlll  2 

T«OPla|.0*Fl 

**■10 

N®2,M* 

8*00*0 

OClTtal.OO 

0CTA*|tlT08 

SLUPfcao.Uf 

GAK*AatO,0 

0*2.8 

h»n.i 

NU2«N/| 

G®*0.78 

6*82.2 

6RP*0«7A 


Of 


100 


109 


C 

C  h I uh  rwco  cutoffpj.o  sec 

C  SPAdAL  ALIASING  CUTOFF*) .8  8CC 

C  NLOwaLU*  FREUUENCV  CUTOFF  NUH«M 

c  NTFN*HICH  FREQUENCY  CUTOFF  NU*i»ER 

C  n|*lFR«8FacUL  aliasing  frequency  cutoff  NUMGF* 

C  FREQUENCY  CUTOFF  NUMBERMlME  SERIES  LCNCTM/CUTQFF  PERIOD 
nLUh*9o 

NYFR*  U9 

NSALFRa)01 

nshi*nsalfr*i 

110  CUNTINUE 


c 

C  initializing  VALUES 
SLFRtQaO.O 
SUDD*0 .0 

110  SHFRtOaO.O 

SUN|N«d,o 
80*1*0, 

SUHJ«o. 

SUNFl-u, 

119  8UNF 2*0  a 

*SUMl*0«0 

*»uh2*0#q 

AVttl*0.0 

AVU2»0.0 


120  PLPOSao.O 

PLNCGao.O 
PLNEIao.O 
8*02*0,0 
UU  20  !■ 1 • N 

129  FiKDan.o 

F2I( ))a0.0 

80  CONTINUE 

00  )u  I* 1 • N02 
FHUDSQ(1)*0.0 

1)0  90  CONTINUE 


Iff 


\  THJS  PORTION  OF  PROGRAM  READS  IN  NAVE  PRESSURE  VALUES  INTO  F|R«F2R  ARRAYS 

C  AND  ASSURES  MATCHING  DAl|  GROUPS  FUR  DlRlCTIUNAL  NAVE  ANALYSIS  OF  TNO 

C  SAGES. 

CALL  8UF  (HSAGEl  .NORTH!  v  MOAT  IfHTIHElt  FI*  •  IDATEI  tlNC1) 

IF(ENO)  GO  TO  1 

CALL  ftuF<"GAGE2«"ON?H*fH0Arf«H? JHt*9FlR  •IDATlIfEND) 

IF ( E NO)  60  TO  1 

IF(lDATS!.EQ.tDA|||)  GO  TO  120 


Figure  2.  Listing  of  main  program. — Continued 
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T 


MO  BACKSPACE  0 

uu  TU  no 

MO  CONTINUE 

IP(HU*6e|.NC,)ii)  CALL  l«ITCH(HCAOtl |R  iFfR  ) 

1«)  426  FUKHATC//*!  GAUGE  NO,  {*#X t  (MONTH  (*  71  *  IDAY  l*  ITIHg  |) 

NRlHU*ll  JMGACf  i*MONTM|#NOATl#H?|Mfl 

■"NlHUiinHGAGC2*N0NTM*fMDAY*tMTlM|2 

11  PUMMAT(I7tl(*X»IT)) 


ISO 


111 


160 


161 


c 

c  THIS  PORTION  or  PROGRAM  CALCULATES  naTER  DEPT M  AT  HAVE  OaCES  AS  mEU  A» 
C  A V t R AGE 8  and  SUM  OF  SQUARES  UP  T]Hf  IfKUl 

00  ui  M1«N 

AVGMAvGl*FtR(  I) 

UZ  AVG26AVG?4P?H(I) 

AVGM*VGl /FLOAT  (N) 

AVG^M AVG2 /FLOAT (N) 

OEPTM«<  AVGMAVG2)/2.tS 

CALL  HfC(OERTN#S«OfeLTT,N,N|ALFR) 

OU  41  I  ■ 1 t  N 
F1M(1)«F|R(  D.AVttl 
F2M(I)>F?R(I)«AVG? 

SUMMSummf  10(1)6*1, 

SUM?*8UM24^2R(  n**2# 

41  L0NT1NUF 

StjHMSuMi/FLUATfN) 

SUHp«suM2/FL0AT  tN) 


170 


c 

c  T*iS  PORTION  OF  PROGRAM  APPLIES  DATA  RlNUU*  TO  TIMt  StRlES«*UATA  *INUU* 

C  VALUES  ARP  REPRESENTED  ttV  Ml) 

00  09  li| #N 

"( 1) "0.9*1 | •0«CHS(T^UP1*PLUAT( 1)/FLUAT(MJ ) ) 


MK(Ui(MM  1)  )♦"<!) 

F2R( l)«(F?R( I)  >**(I) 

*9  CUNT1NUI 


171 


too 


Ml 


C  THIS  PORTION  (IF  PROGRAM  COnPUTtS  SUM  OF  SQUARES  OP  OAtA  nINOUN  MOOIFIEO 

t  1  I ME  SERIES  AS  "ILL  Al  RaTIO  OF  PRE  MJNDUM10  ENERGY  TO  nINDOnED  ENERGY 

00  «1  M1«N 
-SUHMR|UHMF1R<I)  • 

NSUHiaRSUMptFpRd  )**2, 

41  CONTINUE 

*$UHM«SUN1/FL0AT(n) 

^Sun2»«*IUm2/FL0AT  (N) 

RAT IUfaSUM| /«SUM| 

HAT !0?aS0M2/*|u*2 
CALL  FFT (P | Rf P l I tRfO) 

CALL  PPT(P2R«F2I«*»0) 

MEANMFlR(l) 

MEAN2aF2R(l) 


MO 


Ml 


200 


c  This  PORTION  of  PROGRAM  CALCULATES  CO  ANO  QUAD  SPECTRA  VALUES*  AS  "ELL  *1 

C  RAVE  ANGLE  TO  SHORELINE  aNU  ENERGY  CONTRIBUTIONS  OF  EACH  FREQUENCY, 

C  BREAMING  MV(  HEIGHT  AND  BREAKING  NAVE  CELERITY  ARC  ALSO 

C  CALCULATED  IN  THIS  SECTION 

Ml 

00  97  Ja2*N 
F  !H(  J)aP  ftR(if) 

PllU)aPlX(J) 

F2RU)aF|R(J) 

F2imaF2I(J) 

!■  Ml 

97  CONTINUE 
00  96  I a  1 • H 

SUMP  laSUHFMFiR(n**S,tF|l(U9*2, 


Figure  2.  Listing  of  main  program, — Continued 
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2o5 


*«  CONTINUE 

6u*F2munf2*hean2M2. 

IIO  1*0  FU*HAT(//.7X«  Cl(«IOX*  (SlG"A(l)  (,tU«  lENO|Q(IM) 

UU  9*  l«t«*02 

Cti(l)B* lN(t)«FlR( l)tr I  1(1)9^11(2) 
Ut2(I)aFiff(I)9F>l(2)»F2R(l)9F|l(l) 
SlGHA(I)artUAT(I)9TMu*X/(fL0AT(N)9DEk.Tt) 

215  T*1wuP|/ilaHA(l) 

CALL  «VLEN(OiPTH*T*KAH) 

AAaXAH/DtPTN 

1MC12(|).LE. 0.000000001  )  GO  TU  9; 

**On(l  A/(X«t8£  )«AfANrtfl2U2/CI2(Ij) 

220  lf(A»8(PD)#Gr,l eu)  GO  TO 

TMfeTA(n«-ACO|(PO)ABETA 
GO  TU  92 
95  TMtT*m«0,0 
GO  TU  92 

225  93  TffcTA(I)»o,00001 

f «UOSU( X)«F !*(1)**2.+F 1 1(1)992. 

AAb»«KH9H/0EPTH 

J(HK«CUiH(XKH)/CO*H(KAHj 

MU0»UU)aPH0080(I)/(«KP9  92.) 

230  t«GD»OU)aFHQDSOU)*KAUoi 

•0009 SOOOaF  MCDSO ( 1 ) 
a*lTfc(9tl05)l»SlGMA(l)»FHUDIU(l) 

105  FOhmat  (3X' !5t5x  »M2 .9' 7X'Mt«») 

92  CUNT  X Nyt 

235  A*UDiU(  t)9PU(( I )992t9Pl  1(1)992* 

AP9SAKN9B/0CP7H 

«PP9GutH(XAM)/C05H(XKH) 

PPUUtUC I)aFHOO$Q(l)/(lKP**2«) 

F«U0ig(X)apH0O80(I)*PATX0J 
290  MNao,5«<l .♦2#t«KM/|lNMC2t9XKH)) 

CGC  l)a||tGHA(2)90iPTH9PN/XKH 
C(1)«CG(I)/PN 

HG|a(CG(X)92.0aFHODlU(X)9C0t(THETA(l))) 

C  |MG2a0Nf HOPE  CNEPGv  FLU* 

195  )^X«SH01«H62 

•UPCNa|UH|NaPHODSU( X) 

IF(I.CE.NiALFi)  GO  TU  79 
GO  TU  7* 

79  •HMtOatHFRCUAFMODtQU) 

250  76  CONTINUE 

IMI.Lt.NLO*)  GO  TO  77 
GU  TU  7* 

77  •LFREQMLFREG^MUDIQU) 

79  CUNTlNUE 

235  IMI.CE.NVFR)  oo  TO  999 

99  CONTINUE 
999  CONTINUE 

lHC?a|HQ292, 

*N2Tt (4*351 ) 

260  151  FONnaT (  //iftXt  IKMIX'  ItlGNA(I)  (MIX'  (PCT(M4X.  (THETA (  I )  () 

UU  4*  (a) « ND| 

IF ( I *GE .N8ALM)  GO  TO  4* 

PCT«FHOD*0(X)/iUNEN 
IFtPCT .GE.0.025)  GO  TO  99 
265  GO  TU  4* 

49  *MIU(4,*0)If6IQHA(I),PCT#THETA(I) 

30  FONNAT(jXfX5f3(3XtFl6.6)) 

46  CONTINUE 
44  CONTINUI 


Figure  2.  Listing  of  main  program. — Continued 
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*T0 

Cb0<Q9HB/GBP)M,S 
*7  CONTlNut 


*?* 


290 


205 


200 


291 


100 


10* 


110 


11* 


120 


1$ 


1HJS  PORTION  or  PROORAH  MODIFIES  WAVE  GAGE  ANGLE*  TO  BREAKING  WAVE  ANGLES 

ano  computes  longshore  energy  flux  factor* 

DU  9}  UlfNDE 

Jf ( 1 aG£ ,N| aLFR)  GO  TO  090 

)«01n(ThETA( i))*C*/C(I) 
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Figure  2.  Listing  of  main  program. — Continued 
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(1)  Input  data  for  this  program  are  in  the  form  of  digital 
magnetic-tape  records  of  4,100  values,  Vne  first  4  values  of  the 
records  are  the  gage  number,  month,  day,  and  time  of  the  observa¬ 
tions;  the  last  4,096  values  are  the  time-series  pressure  values  of 
the  wave  gage.  In  the  present  program  the  wave  gage  pressures  are 
stored  in  thousandths  of  a  foot  (head)  water  at  0.25-second  inter¬ 
vals.  Subroutine  BUF  reads  time-series  data  into  array  CNTL, 
where  it  is  averaged  to  provide  1,024  time-series  values  of  At  *  1 
second  spacing.  Units  are  also  divided  by  1,000  to  convert  values  to 
feet  (head)  of  water. 

(2)  The  date  groups  of  record  1  and  record  2  are  compared  to 
ensure  that  times  of  records  are  simultaneous;  if  the  times  are  not, 
the  program  searches  the  record  file  until  this  condition  is  met. 
The  two  records  are  than  checked  for  proper  sequence  to  ensure  that 
gage  1  is  analyzed  first.  Subroutine  SWITCH  switches  arrays  if  they 
are  not  in  proper  order. 

(3)  Each  of  the  two  1,024  value  time  series  is  then  analyzed  for 
average  values  which  are  printed  out  along  with  the  average  depth  of 
water  at  each  gage  site.  The  average  value  of  each  of  the  time- 
series  records  is  again  averaged  and  is  added  to  the  height  of  the 
gages  above  the  bottom  to  obtain  the  water  depth: 

DEPTH  «  AVERAGE  1  +  AVERAGES  + 

in  which  AVERAGE  1  is  the  average  of  time  series  1  *  a^O),  AVERAGE  2 
the  average  of  time  series  2  *  a2(0) ,  and  B  the  height  of  sensors 
above  the  bottom. 

An  option  to  apply  a  weighting  function  w( j )  (*  W(I)  in  program) 
has  been  incorporated  before  the  FFT  subroutine  is  called.  In  this 
particular  program  a  cosine  bell  weighting  function  has  been  incorpo¬ 
rated.  If  the  data  window  option  is  selected,  the  two  time-series 
data  records,  which  are  read  into  FIR  and  F2R  arrays,  are  multi¬ 
plied  by  the  following  weighting  function  (cosine  bell) 

W(j)  -  j[i  -  cos 

where  j  is  the  time  step  number  and  N  the  number  of  data  points 
in  series.  If  no  weighting  function  is  desired  in  analysis  set 
w(j)  «  1.0,  which  is  the  Hbox  car”  weighting  function. 

As  the  cosine  bell  function  reduces  the  total  energy  content  of 
the  waves,  the  final  energy  obtained  from  the  FFT  must  be  rescaled  up 
to  the  proper  value.  This  is  accomplished  by  scaling  up  the  time- 
series  pressure  values  by  the  ratio 


<p2> 

<p,2> 


Unwindowed  energy 
Windowed  energy 


■ 


as  discussed  In  equation  (A). 
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(4)  Cospectra  and  quad-spectra  of  the  gages  are  computed  using 
the  following  relationships  (note  in  computer  program  Index,  I  is 
used  for  frequency  counter,  n) : 

Cospectra  «  C12(I)  -  F1R( I)*F2R( I)  +  FI I( I)*F2I( I) 

Quad-spectra  -  Q12(I)  -  F1R( I)*F2I( I)  -  F2R( 1)*F1I( I) 

in  which  FIR  and  FI  I  are  the  real  and  imaginary  parts  of  complex 
transforms  of  time  series  1:  F2R  and  F21  are  the  real  and  imagi¬ 
nary  parts  of  complex  transforms  of  time  series  2. 

(5)  Wave  angle  is  calculated  in  accordance  with  equation  (19). 

0(n)  «  0  -  - ^—5 —  f— i—  .  arctan 

v  '  arcosine  I  k(n) l  C12(n)J 

where  k(n)  is  the  wave  number  calculated  via  linear  wave  theory,  i 
the  spacing  of  gages,  and  6  the  difference  in  alinement  of  gages 
and  shoreline  in  Figure  1. 


Due  to  energy  leakage  problems  in  spectra,  impossible  wave  angles 
can  result  [wave  angles  with  (l/k(n)t  arctan  Q12(n)/C12(n))  greater 
than  1.01.  When  this  happens,  energy  is  lumped  into  a  separate 
category  for  later  scaling  up  of  the  longshore  energy  flux. 


(6)  The  high  frequency  cutoff  in  this  particular  program  has  been 
set  at  2.09  radians  per  second,  which  corresponds  to  a  period  of  3 
seconds  or  NYFR  »  342.  This  value  can  be  reset  in  the  main  program 
by  adjustment  of  NYFR  where 


NYFR 


and  N  is  the  number  of  data  points  in  time  series ,  At  the  spacing 
in  time  of  data  points,  and  Tftp  the  high  frequency  cutoff  period. 
The  spatial  aliasing  frequency  is  computed  in  subroutine  HFC. 


Energy  between  the  spatial  aliasing  frequency  and  the  high  fre¬ 
quency  cutoff  is  put  into  a  special  category  and  used  to  scale  up  the 
final  longshore  energy  flux. 


(7)  Each  frequency  contribution  to  the  onshore  energy  flux  is 
calculated  for  the  gage  site  location  as  follows: 

Onshore  energy  flux  -  2y  | F^Cn ) | 2  Cg(n)  cos  [0(n)] 
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where 


IVn>l 

Cg(n) 

0(n) 

Y 


modulus  of  the  complex  amplitude  spectra  of  wave 
elevation  above  mean  surface  at  gage  site 

group  wave  speed  at  gage  site 

0  *  angle  ->f  wave  direction  (see  Fig.  1) 

specific  weight  of  seawater 


The  onshore  energy  flux  is  then  summed  to  obtain  the  total  onshore 
energy  flux.  In  the  program,  onshore  energy  flux/y  =  HG2. 

(8)  Breaking  wave  height  at  the  shoreline  is  determined  from  the 
mean  square  onshore  energy  flux  via  a  linear  theory  wave  transforma¬ 
tion  formula  which  can  be  simplified  to 

[N/  2  1  /  GB\^  *  ^ 

Ji  16|Fn(n)|2  Cg(n)  cos  0(n)IO.4/_j 


where  GB  is  the  wave  height-to-water  depth  ratio  at  breaking  and  g 
the  acceleration  of  gravity. 

The  choice  of  GB  is  up  to  the  user  although  a  value  of  GB  = 
1.42  has  been  found  by  Komar  and  Gaughan  (1972)  to  best  fit  wave  tank 
data  for  breaking  wave  heights  for  monochromatic  waves.  In  the  pres¬ 
ent  program,  GB  has  been  set  equal  to  0.78  but  can  be  readily 
changed . 

The  breaking  wave  water  depth  is  calculated  from  the  equation 


Hk 


where  d^  is  the  wave  breaking  water  depth  and  GBP  the  ratio  of 
wave  height  to  water  depth  at  breaking. 

In  this  case  a  different  value  of  the  ratio  of  breaking  wave 
height  to  water  depth  can  be  used  in  the  program  for  obtaining  the 
proper  water  depth.  An  assumed  value  of  GBP  *  0.78  from  the  solitary 
wave  theory  in  the  SPM  is  used. 

Linear  wave  celerity  is  assumed  and  breaking  wave  celerity  is 
estimated  as 


The  breaking  wave  height  and  celerity  calculated  in  this  approach 
apply  to  all  frequencies. 
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r  .'If  - 


(9)  Frequency-by-frequency  modif ication  of  wave  angles  is  made 
assuming  linear  wave  theory,  Snell's  law,  and  parallel  bottom  con¬ 
tours  offshore.  The  breaking  wave  angle,  0^(n) ,  is  calculated  from 


®h(n> 


arcsin 


Cb(n)  sin  ©r(n) 

_  ^  _ 


where  the  subscript  r  refers  to  the  reference  gage  location. 

(10)  Longshore  energy  flux  is  calculated  for  each  frequency 
component  (except  the  special  cases  discussed  in  Sec.  II)  using  the 
equation 


p*s(n)  =  Y|Fn(n)|2  Cgb(n)  sin  20b(n) 
and  is  summed  up  to  obtain  a  net  longshore  energy  flux. 


(11)  The  value  of  the  net  longshore  energy  flux  is  multiplied  by 
a  factor  R  which  scales  up  the  total  energy  in  the  spectrum  (below 
the  high  frequency  cutoff).  The  equation  for  scaling  factor  R  is 


R  “  (1  -RTOT) 


where  RTOT  *  RSODD  +  RSHFRQ  when  RSODD  is  the  percent  of  energy  in 
low  frequency  bands  for  which  impossible  values  of  the  cosine  func¬ 
tion  are  calculated,  and  RSHFRQ  is  the  percent  of  energy  between 
spacial  aliasing  frequency  and  high  frequency  cutoff. 

The  final  result  of  analysis  of  the  tw  gage  records  for  trie 
net  longshore  energy  flux  PLNET  is  printed  out,  as  well  as  specific 
frequencies  for  which  impossible  directional  results  occur  and  fre¬ 
quencies  at  which  more  than  2.5  percent  of  the  total  wave  energy  is 
found. 


IV.  SUBROUTINE  DOCUMENTATION 


FFT  Subroutine. 


The  sampled  time  function,  f(j)»  will  be  expressed  as 
N-i 

f(j)  -  l  F(n)  exp(ina>1  jAt) 

n-0 
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in  which 


2n  2it  2tt 

a)  «  -  *  —  *  - 

1  record  length  T  NAt 

tj  *  jAt  *  a  discrete  time  where  j  is  the  integer  time  step 
F(n)  *  a(n)  -  ib(n) 

‘(rh)--*(rh) 

a(0)  =  mean  of  sampled  record 

b(0)  -  b(|)*  0 

Because  negative  indexes  are  not  readily  handled  by  most  FORTRAN  compilers, 
the  summation  extends  over  the  interval  0  <  n  <  N  -  1,  rather  than  over  the 
symmetric  interval  -  N/2  <  n  <  N/2.  From  the  definition  of  the  coefficients 
above,  it  is  clear  that  the  coefficients  a(n)  and  b(n)  for  n  >  N/2  contain 
no  additional  information . 

The  inverse  relationship  completing  the  FFT  pair  is 

i  N 

F(n)  «—  l  f(j)  exp(-inu)  jAt) 

N  J-l 

As  an  enumeration  of  the  complex  FFT  coefficients,  suppose  the  series  of  8 
values  is  considered,  N  =  8.  The  coefficients  would  be 

F(0)  =  a(0) 

F(l)  -  a( 1)  -  ib(l),  F(7)  -  a(7)  -  ib(7)  -  a(l)  +  ib(l) 

F( 2)  -  a( 2)  -  ib( 2) ,  F(6)  -  a(6)  -  ib(6)  -  a(2)  +  ib(2) 

F(3)  -  a(3)  -  ib(3) ,  F(5)  -  a(5)  -  ib(5)  -  a(3)  +  ib(3) 

F(4)  -  a(4) 

This  pattern  prevails  for  all  sets  of  FFT  coefficients,  regardless  of  the 
value  of  N.  Both  F(0)  and  F(N/2)  are  real  and,  as  noted  previously,  the 
coefficients  F(n)  for  n  >  N/2  really  contain  no  additional  information. 
The  FFT  subroutine  used  here  requires  that  the  number  of  data  points,  N, 
provided  be  an  integral  power  of  2,  i.e., 

N  -  2k 

where  K  is  an  Integer.  Thus  analyses  of  512,  1,024,  or  2,048  data  points 
(K  -  9,  10,  11)  would  be  suitable  with  this  subroutine. 


L _ . 


23 


The  fol Lowing  two  requirements  are  satisfied  in  the  FFT  subroutine. 


(a)  By  operating  on  the  sampLed  function,  obtaining  the  F(n) 
coefficients  and  carrying  out  the  inverse  FFT  (FFT-1),  the  original 
time  function  is  recovered.  Schematically, 


f(j)  - 


(b)  The  mean  square  of  the  sampled  time  function  is  equal  to  the 
sum  of  the  squares  of  the  moduli  of  the  FFT  coefficients,  F(n) ,  i.e., 

2  N  N-l 

-  I  -  l  |F(n)|2 

N  j=l  n=o 

a.  Calling  Statement:  SUBROUTINE  FFT  (FR,  FI,  K,  ICO)  (see  Fig.  3).  FR, 
FI  »  real  and  imaginary  coefficients  in 


FFT 


F(n) 


FFT 


-1 


f(j) 


F(n)  -  FR(n)  -  iFI(n) 

K  =  power  of  two  (i.e.,  N  *  2K) 


ICO  =  control  whether  FFT  or  (FFT)-1 


operation  is  desired  if 


ICO 


0  FFT 
1  -►  (FFT) 


-1 


When  entering  the  subroutine,  FR  is  the  time  sequence  f(j)  and  FI  is 
arbitrary.  When  exiting  the  subroutine,  FR  and  FI  are  the  real  and  imagi¬ 
nary  parts  of  the  complex  transform,  respectively;  e.g.,  input  is 

K  -  5 

ICO  -  0 

f(j)  -  l.o  +  2.0  cos  +  3.0  cos 


-  0.6  sin  2^At)  -  1.4  sin 
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Figure  3.  Listing  of  FFT  subroutine* 


Data  Input  to 

Program . 

f(j)  values 

at  6.0'JO 

5.080 

3.750 

2.184 

At  *  1  second 

(’ .  590 

-0.829 

-1.900 

-2.506 

intervals 

-2.600 

-2.215 

-1.451 

-0.465 

(32  values) 

0.562 

1.445 

2.034 

2.229 

2.000 

1.391 

0.513 

-0.475 

-1.390 

-2.054 

-2.322 

-2.109 

-1.400 

-0.257 

1.188 

2.755 

4 .238 

5.438 

6.189 

6.386 

FR  » 

6.000 

5.080 

3.750 

2.184 

(32  values) 

0.590 

-0.829 

-1.900 

-2.506 

-2.600 

-2.215 

-1.451 

-0.465 

0.562 

1.445 

2.034 

2.229 

2.000 

1.391 

0.513 

-0.475 

-1.390 

-2.054 

-2.322 

-2.109 

-1.400 

-0.257 

1.188 

2.755 

4.238 

5.438 

6.189 

6.386 

FI  - 

0.000 

0.000 

0.000 

0.000 

(32  values) 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

Calling  Statement : 

FFT  (XR, 

XI,  5,  0). 

Output : 

1.000 

1.000 

1.500 

0.000 

a(n)  coefficients 

0.000 

0.000 

0.000 

-0.000 

(32  values) 

-0.000 

-0.000 

-0.000 

-0.000 

-0.000 

-0.000 

-0.000 

-0 .000 

-0.000 

-0.000 

-0.000 

-0.000 

-0.000 

-0.000 

-0.000 

-0.000 

-0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

1.500 

1.000 

b(n)  coefficients 

0.000 

-0.300 

-0.700 

-0.000 

(32  values) 

-0.000 

-0.000 

-0.000 

-0.000 

-0.000 

-0.000 

-0.000 

-0.000 

-(K900 

-0.000 

-0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.700 

0.300 

At  (time  step)  *  1 

second  in 

above  example 

. 

2.  HFC  Subroutine, 

rhis  subroutine  resets  the  spatial  aliasing  frequency  cutoff  to  a  higher 
frequency  than  would  be  the  case  for  normal  incidence  of  waves  to  gage  pair. 
In  the  present  version  of  this  subroutine,  it  has  been  assumed  that  the  maxi¬ 
mum  angle  which  the  wave  crests  can  make  with  the  gage  pair  axis  is  45°.  The 
spatial  aliasing  criteria  are  expressed  in  Figure  1,  where  for  proper  resolu¬ 
tion  of  wave  direction  the  following  criteria  must  be  met 

fccos  [Q(n)  -  el  <  i 
k(n)  Acos  (0(n)  -  01  <  k(n) 
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The  proper  spatial  aliasing  frequency  to  correspond  with  the  special  aliasing 
wave  number  cutoff  is  found  from  the  normal  wave  dispersion  relationship. 


Calling  Statement:  HFC  (DEPTH,  S,  DELTT ,  N,  NSALFR)  (see  Fig.  4). 

DEPTH  “  depth  of  water  at  gage  site  (from  maJ-^  program) 

S  *  spacing  of  wave  gage  pair  (from  main  program)  (■  l  in  text) 

DELTT  *  time-step  increment  between  values  in  time  series  analyzed 
(from  main  program) 

N  *  exponent  of  2  describing  number  of  time  series  values 

(from  main  program) 

NSALFR  *  integer  number  for  spatial  aliasing  frequency  cutoff 


t 


s 


10 


c 

c  SUBMUUTlNfc  nrClHlGH  FREUUtNCY  CUT  OF  F /SP  At  I  *L  ALIASING  FRtOUENCV) 
C  RESETS  ALIASING  CUTOFF  To  HIGHER  FHEOUlNCY 

C  am 0  UN  ASSURED  HAXJHUM  NAYE  ANCLE 

SUttRUUTlNI  HFC(OE^TH<»fOELTT»NtNSAL P«> 

C  8PAC1AL  ALIASING  ASSURES  NaVE  ANCLES  LESS  THAN  #1  016»UI 
XAS*J,l«iS9/0.707 
XA Ha AKSA DEPTH/8 

ai6SUa)|,|*(XKH/DEPTH)ATANH(XKH) 

9 1 CHE* SORT (1 1080) 

HECLN«FLOAT(N)ADILTT 

NSALFHaSXCHFARECLN/S.teJ 

RETURN 

(NO 


Figure  4.  Listing  of  HFC  subroutine. 


3.  SWITCH  Subroutine. 


This  subroutine  is  set  up  to  interchange  time-series  data  arrays  in  the 
instance  when  gage  2  data  are  processed  before  gage  1  data  (see  Fig.  5).  If 
the  first  gage  record  processed  is  not  equal  to  the  appropriate  number  of  the 
gage,  as  specified  in  main  program,  data  arrays  of  first  and  second  gage 
records  are  interchanged. 


i 


s 


10 


18 


C 

C  6U8NUIITINE  SWITCH  EXCHANGES  LOCATIONS  of  TINE  SERIES  DATA  TO  ASSURE 
C  0AGC1  IS  STORED  IN  FIRST  ARRAY  ANO  CAGES  IN  SECOND 

SUBROUTINE  Srt!TCH(Ml*M2'FlRiF2K) 

0 1  HENS  ION  F1R(10I«)'F2R(102«)'F3R(|02R) 

R1»mi 

Ht«HJ 

H2RMJ 

00  10  I«1'102« 

FJN(1J,F1R(I) 

F1H(I)rF|R(I) 

F2H<I)«F3R(I) 

10  CONTINUE 
RETURN 
I  NO 


Figure  5.  Listing  of  SWITCH  subroutine. 


4.  WVLEN  Subroutine. 


This  subroutine  accepts  wave  period  and  water  depth  as  input  and  calcu¬ 
lates  wave  number  as  output  via  a  Newton-Raphson  iteration. 
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Calling  Statement:  WVLEN  (DPT,  PER,  XKH)  (see  Fig.  6). 

DPT  »  water  depth  (from  main  program) 

PER  *  wave  period  (from  main  program) 

XKH  *  wave  number  *  water  depth  (calculated  in  subroutine) 


t 


$ 


to 


if 


10 


C  LENGTH  ITERATION  |U*ROUTXNI.,thj*  buyout  ini  CALCULATE*  NAVILENGTh 

C  VIA  n|*70N*paPHIon  ITERATION  Ul JNO  RERI06*NATER  DEPTH  INPUT 

C  REHmwm  Rl«XOO 

C  0Pt«HATE»  depth 

C  XftHatAVf  a ATER  DEPTH 

•UHRUUTINC  »*VUN(DPTtPER»XRH) 
kRHu«(6.t83t«9)/PER)PP2PDPT/lt»a 
2P(XHH0-A.J)a»l ti 
\  XKHBXKHO 
GO  TO  f 

2  XKHa9UHT(XKH0) 

I  SHa8JNH(XNH) 

CHpCUtHCNKH) 

EP*aXAH0»XHHA$H/CH 

KbPtBaXHH/CHMtalHSC* 

UV*H«.tP*/*LQPE 
IP (Att*(DXKH/XKH)«O90Q0l )*»*•* 

«  XKhaHKHADIRH 
GO  TO  1 
9  CONTINUE 
RETURN' 

END 


Figure  6.  Listing  of  WVLEN  subroutine. 


5.  BUF  Subroutine. 


This  subroutine  is  set  up  to  read  in  wave  gage  files  from  magnetic  tape. 
The  data  records  consist  of  arrays  of  4,100  values,  the  first  four  of  which 
are  the  gage  number,  month,  day,  and  time  of  wave  record.  The  remaining  4,096 
values  represent  pressures  in  thousandths  of  a  foot  (head)  water.  The  data 
are  returned  to  main  program  as  a  wave  gage  number-date  series  and  a  time 
series  of  4,096  values  of  pressure  in  feet  (head)  of  water.  Two  records  are 
processed  in  one  pass. 

Calling  Statement:  BUF  (MGAGE,  MONTH,  MDAY,  MTIME,  CNTL,  1DATE ,  END)  (see 
Fig.  7). 

MGAGE  -  number  of  gage  (read  from  tape) 

MONTH  ■  month  of  observation  (read  from  tape) 

MDAY  -  day 
MTIME  *  time 

CNTL  *  control  array  of  4,096  pressure  values  in  feet  (head)  of  water 
returned  to  main  program 

IDATE  *  summed  time  group  for  time  comparison  between  gages 
END  *  logical  end 
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1 


5 


10 


15 


io 


2* 


10 


1« 


40 


C 

t  SUBROUTINE  HUF  READS  In  nAVt  GAGE  DATE  INFU  AND  TIHE  •CRIES  DYNAMIC 

C  RRESSURf  VALUfl  IN  FECI  HEAD  Of  hATER 

C  THIS  SUBROUTINE  READS  4004  TIME  SCRIES  VALUES  AND  AVERAGES  TO  OBTAIN 

C  1024  VALUES  FOR  MAIN  PROGRAM  ANALYSIS 

C  "GAGfciGAGF  NUMRER 

C  munThhhunTh  UF  RECORDING 

C  MQAVBftAV  OF  RECORDING 

C  MtlMfeiTIHI  OF  RECORDING 

C  HE AL* ARRAY  OF  AVERAGED  TIME  SERIES  VALUES 

SUttHUUTINF  SUF ( MGAOE *  MONTH, M^ATtMTIMEtRCALtlOATCtlND) 

DIMENSION  CNTL(40S4)f U(SOOO) 

DIMENSION  REALd024) 

LOGICAL  END 
DO  )2  J* 1 , 4044 
CNTLIJ)«0.0 
12  CONTINUE 

BUFFER  lNfG,l)(!A(l)»lA(S000I) 

IF(UNIT(4) )l0,20t)0 

20  r«int  ii,nA(n,i«i,4i 

11  FUHMATn  RARITY  ERROJ  ON(tAlT) 

10  CONTINUE 
HGAGtilA(l) 
month* 1  a ( 2) 

MOAY4IACJ) 

MTIHt«IA(4J 

IOATfc*IA(2)M4Cm!AU) 

M4 

UO  2)  J* 1 , 4044 
MM  | 

CNTL(J)«IA(K) 

25  CNTLU)iCNTL(J)/1000, 

00  24  j440SS*4044 

24  CNTLU)4CNTL(40S7) 

J«l 

00  27  L«1#I024 

RCALiL)4CCNTL(J>ACNTl<JM)4CNTL(j4|)ACNTL(JAl>)/4. 

J«J*A 

17  CONTINUE 
RETURN 

10  End* • TRUE • 

RETURN 

End 


Figure  7.  Listing  of  BUF  subroutine* 
V.  SAMPLE  OUTPUT 


Three  examples  of  output  are  presented  for  different  dates  for  the  wave 
gage  pair  at  Channel  Islands  Harbor  (Fig.  8).  The  year  the  data  was  taken  was 


1975. 


The  first  set  of  frequencies  lists  amplitude  modules  squared  of  wave  data 
having  impossible  direction  results.  The  sum  total  of  this  energy  (in  decimal 
percent)  is  listed  as  the  quantity  RSODD  in  the  variable  output  at  the  bottom 
of  the  output.  In  the  case  of  the  wave  data  taken  on  7-26-1600,  the  inco¬ 
herent  data  amounted  to  0.004  (0.4  percent)  of  the  total  energy  in  the  wave 
record. 

The  second  set  of  frequencies  listed  provides  the  wave  direction  for  the 
frequency  bands  having  a  significant  part  of  the  energy  (>2.5  percent).  In 
the  case  of  the  wave  record  taken  on  7-26-1600,  it  is  seen  that  the  wave  angle 
is  reasonably  consistent  from  the  frequency-to-frequency  band  and  is  approxi¬ 
mately  0.70  radian  (40.1°). 

The  variable  list  provided  at  the  bottom  of  the  sampled  output  gives 
values  of  most  Importance  in  the  analysis  of  wave  information  for  longshore 
energy  flux.  The  longshore  energy  flux  output  is  in  pounds  per  second  units; 
the  output  in  the  first  example  is  89.23  pounds  per  second. 
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Example  1 


GAUGE  NO. 

HUNTH 

OAT 

T1M9 

Ill 

T 

26 

1609 

ii  a 

T 

26 

1699 

i 

II9"A(l) 

7*080(1) 

i 

•oi**o* 

.000019 

4 

•9)«9«* 

.00001* 

5 

•910689 

.0000(7 

7 

•0**991 

•0000)6 

9 

.09922) 

.00001) 

10 

•081)99 

•0000)1 

11 

•OOT999 

.0000*8 

14 

•08990) 

•  000909 

24 

•1*7262 

•0001)7 

25 

»l9))98 

.0000*1 

27 

•169670 

.09000) 

24 

•1779*2 

.000099 

SO 

.18*078 

.0000*0 

12 

.1*6)90 

.000069 

1) 

.202*89 

.0000)9 

42 

.297709 

•  00001* 

45 

•>76117 

.0000*1 

65 

•>989)9 

•  0000*) 

1 

)16H*(!) 

67 

•*11106*9 

60 

•*172*277 

71 

•  **79 2M9 

7« 

•*9*099)1 

75 

•*601**)* 

70 

•*7860201 

74 

.*8«T)79) 

2CT 

•01404604 

loMIKtl 

•10)9919* 

•048907*0 

.0290*019 

ituirii! 


TMETA(I) 

•7017699) 

•rr*M«M 

•MTmi) 

•69199*90 

•7|«|))9« 

•umiH 


N8*C7R9 

20! 

0)77*  07 

*AtE»  AT  GAUGE  0|Te» 

2S.I 

*V8|* 

#1.411 

AVG24 

14,444 

9U*i* 

•  224 

0U*2« 

•  214 

M9UH1* 

.004 

N0UH0P 

•  006 

RATIO]* 

2*710 

HATIQ24 

2,724 

IU*E** 

•10411410 

ORCaKINO 

»Ayt  HEIGHT  H04 

1*01 

8R£aK JR6 

HAyE  CEliGITy  C04 

11,10 

R6000* 

•  0040 

R)h7MQ* 

•  2411 

MLMtfs 

•  0179 

7L709* 

44*6021 

7LNE8* 

•§•4150 

7U*ET* 

04*2271 

Figure  8. 

Three  examples  of  output  for  wave 

gage  pair  at 

Channel 

Islands  Harbor 

• 

Example  2 


Figure 


gauge  no 

.  MONTH  OAT 

TJHI 

Ml 

T  26 

1600 

312 

T  2* 

1800 

1 

211***111 

7*060(1) 

1 

.00*1 J* 

•000160 

2 

•012272 

.000099 

3 

•010*08 

•000007 

5 

•010*00 

•  000013 

6 

,0**007 

•  000066 

9 

•005221 

•  000164 

10 

,06115* 

•000036 

1  1 

•067**5 

•  oooooo 

13 

•07*767 

•  000166 

1  4 

•005*01 

•  000201 

1? 

.0*201* 

•  000137 

1* 

.0*0175 

•  0001  14 

10 

.110**7 

•  000055 

1  4 

•116501 

•  000061 

21 

.1*1126 

•  000004 

26 

.15*51* 

•000072 

26 

•171006 

•000006 

30 

,10*070 

•000020 

31 

•1*021* 

•000006 

«2 

.15770* 

•000026 

56 

,15500* 

•000064 

I 

IICHA(IJ 

PCT 

rNCTA(i) 

64 

.*172*277 

•02575219 

•75201011 

71 

•*601**1* 

•04645516 

•65061227 

76 

•*6611016 

•01141793 

•6*562017 

77 

.*72*660* 

•042)7300 

•71**0726 

NIAt fHm 

203 

DEPTH  op 

M*TtH  AT  GAUGE  one* 

2«.0 

A  VO  |  ■ 

2?.2*6  A V02* 

20,606 

SUM}! 

.29*  SUh2« 

•  29) 

«$UM1« 

.00*  N8UH2* 

,076 

P  A  T  1 0  J  ■ 

1.57*  HAT1U2* 

3*741 

9UH(Nb 

,1101**70 

0HC A* 1NG 

M Ay E  HEIGHT  HH* 

3*61 

8KEi*lNG 

"*vl  cututfr  c0* 

12,22 

96000* 

■00*0  RSHFHOa 

•  3742 

P|UF Rui  • 0 1 1 

PLPQ8* 

125.7115  PI.NE6* 

-25,6734 

PLNtT« 

10ft, 0*01 

8.  Three  examples  of  output  for  wave  gage  pair  at  Channel 
Islands  Harbor, — Continued 
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Example  3 


gauqe  no. 

month 

0  At 

TIME 

Ill 

I 

26 

2000 

na 

7 

26 

2000 

i 

*i«h*;  i) 

FHOIO(I) 

t 

.006114 

.0009)0 

i 

•012272 

,000)47 

j 

,010*00 

.000170 

* 

•024544 

.000052 

y 

.042951 

.00001) 

6 

.0*0007 

•OOOOtl 

10 

.061159 

,00001* 

11 

,ori6ii 

.ooooT* 

11 

,079767 

•000085 

It 

•10*111 

.0000*0 

19 

.11650) 

,000008 

21 

•1*1126 

.000011 

12 

.19*150 

.000009 

1« 

.100*21 

.000100 

IS 

.11*797 

.00010* 

16 

.22009) 

,000019 

41 

.2*7700 

*  ?  0  0  0  4  4 

ss 

.1)7476 

.000116 

71 

.4)5*51 

.001082 

1 

IIOMA(I) 

FCT 

NMlPQa 

202 

oipth  OF 

hateh  at  oa mqi  lire* 

21.0 

AVCj« 

•  1 1 702 

AV92i 

20.287 

lUHi* 

•  291 

»Uh2« 

.246 

nium 1 n 

.071 

*IUM*R 

,081 

RAT  101* 

1.R60 

RATIU26 

1,294 

IUmen« 

•19101917 

ORE aK 1  NO 

«AvE  HEJ9HT  H06 

1.64 

0HEAKINQ 

MAyl  CeiBRlTY  CO* 

12.24 

R6000* 

,,009* 

RlHFHUa 

.9*75 

PtROl* 

1)5.9)67 

RUNES* 

-40.22*7 

runet* 

*5.1070 

8.  three  examples  of 

output 

for  wave  g 

TMfTA(l) 


>0141 


Islands  Harbor. — Continued 
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